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We consider a system of spins on the sites of a three-dimensional pyrochlore lattice of corner-sharing tetra- 
hedra interacting with a predominant effective xy exchange. In particular, we investigate the selection of a 
long-range ordered state with broken discrete symmetry induced by thermal fluctuations near the critical region. 
At the standard mean-field theory (s-MFT) level, in a region of the parameter space of this Hamiltonian that 
we refer to as Rs region, the ordered state possesses an accidental 17(1) degeneracy. In this paper, we show 
that fluctuations beyond s-MFT lift this degeneracy by selecting one of two states (so-called t /)2 and iff) from 
the degenerate manifold, thus exposing a certain form of order-by-disorder (ObD). We analytically explore this 
selection at the microscopic level and close to criticality by elaborating upon and using an extension of the so- 
called TAP method, originally developed by Thouless, Anderson and Palmer to study the effect of fluctuations in 
spin glasses. We also use a single-tetrahedron cluster-mean-field theory (c-MFT) to explore over what minimal 
length scale fluctuations can lift the degeneracy. We find the phase diagrams obtained by these two methods to 
be somewhat different since c-MFT only includes the shortest-range fluctuations. General symmetry arguments 
used to construct a Ginzburg-Landau theory to lowest order in the order parameters predict that a weak magnetic 
moment, rriz, along the local ( 111 ) (z) direction is genetically induced for a system ordering into a ip 2 state, 
but not so for ijjs ordering. Both E-TAP and c-MFT calculations confirm this weak fluctuation-induced 
moment. Using a Ginzburg-Landau theory, we discuss the phenomenology of multiple phase transitions below 
the paramagnetic phase transition and within the F 5 long-range ordered phase. 


I. INTRODUCTION 

In the study of condensed matter systems, mean-field 
theory' is often the simplest starting point to obtain a quali¬ 
tative understanding of the essential physics at play prior to 
carrying out a more sophisticated analysis. A standard mean- 
field theory (s-MFT) replaces the many-body problem with a 
simpler problem of a one-body system interacting with an av¬ 
eraged field produced by the rest of the interacting particles. 
In systems with competing or frustrated interactions, s-MFT 
may yield a number of states with a degenerate minimum free 
energy below the mean-field critical temperature, If 

these degeneracies are accidental, that is not imposed by ex¬ 
act symmetries of the Hamiltonian, they may be lifted by the 
effects of thermal"' or quantum^ fluctuations, a phenomenon 
known as order-by-disorder (ObD)."*"^ 

The concept of ObD was originally proposed by Villain 
and collaborators as the ordered state selection mechanism 
for a two-dimensional frustrated Ising model on the domino 
lattice."' Since this seminal work, ObD has been theoretically 
identified and discussed for many highly frustrated magnetic 
models.^'* These systems generically possess an exponen¬ 
tially (exp[W“]) large number of classical ground-states (here 
N is the number of spins in the system and a < 1). As a re¬ 
sult, highly frustrated magnetic systems are intrinsically very 
sensitive to fluctuations or energetic perturbations. In the con¬ 
text of experimental studies of real materials, it is difficult to 
distinguish a selection of an ordered state via fluctuations from 
one that would arise from energetic perturbations beyond the 
set of interactions considered in a restricted theoretical model. 
Consequently, undisputed examples of ObD in experiments 


have remained scarce.'® 

Quite recently, following an original proposal going back 
ten years,and building on an earlier study, ^"' several 
papershave put forward compelling arguments for ObD 
being responsible for the experimentally observed long- 
range order in the insulating rare-earth pyrochlore oxide^^ 
Er 2 Ti 207 . In this compound, Er''+ is magnetic and Ti^+ 
is not. The key observation in those works is that the ac¬ 
cidentally degenerate classical ground states are related by 
operations with a (7(1) symmetry.'® This set of classically 
degenerate states form the so-called Fs manifold.^^ For a 
range of interaction parameters of the most general symmetry- 
allowed bilinear nearest-neighbor pseudospin 1/2 exchange¬ 
like Hamiltonian (referred to as fi in Eq. (la) below ) on the 
pyrochlore lattice of corner-sharing tetrahedra (see. Eig. 3), 
the (7(1) degeneracy is exact at the s-MET level as long as 
the cubic symmetry of the system remains intact. As a short¬ 
hand, we henceforth refer to this region of exchange parameter 
space as the r 5 region, as referred to in the Abstract. Refer¬ 
ence [16] showed that the (7(1) symmetry is robust against 
a wide variety of perturbations added to H and argued that 
essentially only fluctuations can efficiently lift the degener¬ 
acy in Er 2 Ti 207 when considering bilinear anisotropic inter¬ 
actions of arbitrary range between the pseudospins.Simi¬ 
lar arguments were made in Ref. [15]. In this compound, 
a particular long-range ordered state, the state (see Eig. 
(1)), 22,23,26,28-31 jg rejected. Considering H on the pyrochlore 
lattice, Wong et a/.'^ and Yan et studied the effect of 
quantum'^’^^ and thermaP^ fluctuations and established a gen¬ 
eral phase diagram for this Hamiltonian at T = O®'. 

The investigations reported in Refs. [14-17, and 32] fo- 
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cused on identifying the mechanism of ground state selec¬ 
tion by taking into account the harmonic quantum'^'^’^^ or 
classicalspin fluctuations about a classical long-range or¬ 
dered state. On the other hand, the problem of state selection 
at temperatures near the critical transition temperature to the 
paramagnetic phase has received significantly less attention. 
Although there have been numerical studies of ObD selection 
of the -02 state at T < 21^04,15,22,23,32-34 Qj. ^pQjj approach¬ 
ing Tc from above, an analytical study specifying the role of 
the individual microscopic anisotropic spin-spin interactions 
in the ObD mechanism at T ~ has, to the best of our 
knowledge, not yet been carried out. The problem of selec¬ 
tion at T <Tc is not only of relevance to the phenomenology 
of Er2Ti207 or other pyrochlore magnetic compounds,but 
it is of considerable interest for all highly frustrated magnetic 
systems proposed to display an ObD mechanism. 

There are situations where a sort of ObD occurs near Tc 
which differ from the textbook cases^^’^ where the state se¬ 
lection near T = 0 ^, proceeding either via thermal or quan¬ 
tum selection, is leveraged upon all the way to the long- 
range ordered state selected at Tc- For example, the long- 
range ordered state selected by ObD can in principle be dif¬ 
ferent for the T = 0 + and T < Tc regimes. This occurs, 
for example, for classical Heisenberg spins on the pyrochlore 
lattice interacting via nearest-neighbor antiferromagnetic ex¬ 
change and indirect Dzyaloshinskii-Moriya interaction.In 
another class of problems, different competing long-range or¬ 
dered states may have the same free-energy at the s-MFT level 
only over a finite temperature interval, T* < T < Tc with 
a transition to a non-degenerate classical ground state at T*. 
There, thermal fluctuation corrections to s-MFT can select one 
of the competing states over the T* < T < Tc window. This 
is what is predicted to occur in the multiple-fc state selection 
in the Gd2Ti207 pyrochlore antiferromagnetic between 0.7 K 
and 1.0 Perhaps the most exotic cases arise when a relic 
of ObD occurs at the critical temperature, while the classical 
ground state is not degenerate for the Hamiltonian considered, 
but would be for a closely related Hamiltonian minus some 
degeneracy-lifting weak energetic perturbations.At the 
phenomenological Ginzburg-Fandau level description, these 
cases are not paradoxical. They only become so when one is 
trying to ascribe a microscopic description to the physics at 
stake and the origin for, at least, one further equilibrium ther¬ 
modynamic phase transition at some temperature below the 
paramagnetic transition at Tc- Finally, and generally speak¬ 
ing, one may ask whether a discussion of ObD at T = O'*' 
is of pragmatic usefulness given that most experiments for 
which ObD may pertain typically proceed by cooling a mate¬ 
rial from a paramagnetic disordered phase to an ordered phase 
and not going below a certain baseline nonzero temperature 
constrained by the experimental set-up. This broader context 
provides the motivation for our work which aims to go be¬ 
yond the sole consideration of a phenomenological Ginzburg- 
Fandau theory that contains all terms, relevant or not (in the 
renormalization-group sense), allowed by symmetry and to 
expose how the different competing microscopic interactions 
participate to the degeneracy-lifting near the transition at Tc- 

In the temperature regime near a phase transition, the har¬ 


monic approximation describing low-energy excitations usu¬ 
ally employed in theoretical discussions of ObD^“* at T = 0"^ 
is not physically justified since large fluctuations typically 
accompany the phase transition to the paramagnetic state. 
From a fundamental viewpoint, it is thus highly desirable to 
study the role of fluctuations beyond s-MFT using a different 
method that can be applied at temperatures close to the critical 
region. A possible route to tackle this problem was paved by 
Thouless, Anderson and Palmer (TAP) in their study of the ef¬ 
fect of fluctuations in spin glasses.^® For the case of the Ising 
spin glass model, Plefka showed that the TAP correction can 
be systematically derived using a perturbative expansion;^® an 
approach that was later generalized by Georges and Yedidia 
who calculated higher order terms in the Plefka’s perturbative 
expansion."^' This approach, which we refer to as extended 
TAP (F-TAP), consists of a high-temperature expansion in 
(3 = l/k^T about the s-MFT solution. It captures correc¬ 
tions to the s-MFT free-energy by including fluctuations in 
the form of on-site linear and nonlinear susceptibilities. A 
second approach allowing one to go beyond s-MFT is the clus¬ 
ter mean-field theory (c-MFT)."^^’"^^ This approach treats short- 
range fluctuations within a finite cluster exactly while taking 
into account the inter-cluster interactions in a mean-field fash¬ 
ion. 

In the present work, we focus on the problem of ObD 
selection near criticality in a spin model on the pyrochlore 
lattice with predominant xy interactions using the F-TAP 
method. We have recently used this method to investigate 
the problem of partial multiple-fc order in pyrochlore mag¬ 
nets.^® In this paper, we also employ c-MFT. We concentrate 
on the r5 manifold for the T < T^^ regime since Fs is not 
only an interesting theoretical playground according to recent 
investigations,it is also of potential relevance to 
real pyrochlore materials, such as Fr2Ti207, proposed to dis¬ 
play an ObD mechanism. In addition, the interesting 

case of distinct ObD selection at T = 0 + and Tc, reported in a 
pyrochlore system with anisotropic spin-spin coupling,fur¬ 
ther motivates us in investigating the role of the various spin 
interactions in the state selection at T < in a general model 
of interacting spins on the pyrochlore lattice. 

The rest of the paper is organized as follows. In Section II, 
we present the bilinear nearest-neighbor spin model defined 
by the Hamiltonian T-L of Fq. (la) on the pyrochlore lattice 
and discuss its symmetries. Focusing on the Fs manifold, we 
present a Ginzburg-Fandau (GF) symmetry analysis to spec¬ 
ify the general form of the lowest-order anisotropic terms al¬ 
lowed in the GF free-energy (J^gl) that can lift the accidental 
C/(l) degeneracy found in s-MFT. We show that the fluctua¬ 
tion correction terms to the s-MFT free-energy select either 
the 02 or the 03 long-range ordered state of the F5 mani¬ 
fold (see Fig. 1 and Appendix A for definition) as well as 
induce a weak moment ruz along the local [111] direction for 
the 02 state. In Section III A, we present the F-TAP method 
and determine the phase boundary at Tc between the 02 and 
03 states in the space of anisotropic spin-spin coupling con¬ 
stants. We investigate in Section IIIB how short-ranged fluc¬ 
tuations acting on the length-scale of a single-tetrahedron can 
lead to ObD when incorporated in the self-consistent scheme 
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FIG. 1. (Color Online). Spin configurations of the (a) non-coplanar 
tp2 State, as observed in Er 2 Ti 207 ,^® and the (b) coplanar tps state 
on a single tetrahedron. The pattern is the same on all tetrahedra 
in a pyrochlore lattice which is a face-centered cubic lattice with a 
tetrahedron basis and these two states, which belong to the F 5 man¬ 
ifold, are thus said to have a k — 0 ordering wave vector. The spin 
(in blue) at the a sublattice points along (a) the Xa axis of the local 
[ 111 ] reference frame for the ip 2 state and along (b) the ija axis of 
the same local frame for the t/is state. The local Za axis at sublattice 
a points along the local cubic [ 111 ] axis at that site such that the lo¬ 
cal Xa,ya, Za triad of orthogoual unit vectors (orange arrows) fulfills 
Xa y- Va = Za- As illustrated in panels (a) and (b), the local Za axis 
for sublattice a points directly out of the tetrahedron primitive basis 
cell. See Appendix A for more details. 


four pyrochlore sublattices (see Fig. 1). Si can be treated 
classically as a 3-component vector or quantum mechanically 
as an operator such as a Si = 1/2 pseudospin. In the context 
of magnetic rare-earth pyrochlores, Si would represent either 
the total angular momentum J within a simplified model of 
J — J coupling,or a pseudospin 1/2 describing the 
single-ion ground state doublet. 


As we are foremost interested in the selection of classical 
ordered phases at 0 ^ T < Tc, we shall treat Si generally 
classically with [Sil = 1/2 for all i. However, in Section 
IIIB, where we use the c-MFT method, we consider Si = 
1/2 quantum mechanically, mostly for computational ease. 
In Eq. (1), Jzz, J±, J±± and Jz± are the four symmetry- 
allowed independent nearest-neighbor exchange parameters, 
while Qj = — Jij are bond-dependent phases on a single tetra¬ 
hedron defined in Refs. [16, 45, and 49].^° 


It is important for the discussion that follows to split into 
the two terms T-Lq and TLi in Eq. (la) and consider the sym¬ 
metry properties of each term. T-Lq has a (7(1) symmetry: it 
is invariant under a rotation of Si by an arbitrary global angle 
about the local (111) axes. On the other hand, the TLi term in 
Eq. (la) is only invariant under rotations of Sihy reducing 
the symmetry of "H to Zq (Ca)!!!] x Z2). 


of c-MET. The c-MET method also allows us to explore semi- 
quantitatively the role of nonzero effective Ising exchange 
(Jzz) on the selection at Tc- In Section IV, we briefly discuss 
the possibility of multiple phase transitions in xy pyrochlore 
magnets as temperature, for example, is varied. Einally, we 
close with a discussion in Section V. A number of appendices 
are provided to assist the reader with the technical details of 
the E-TAP calculations. In Appendix A, we detail the spin 
configurations of the 1/12 and -03 states. In Appendix B, we 
analyze the symmetry properties of the ■02 and the i/a states 
to show that a fluctuation-induced local moment is only 
compatible with a 02 long-range ordered state. We provide 
in Appendix C the details of the E-TAP calculations and the 
diagrammatic approach employed. 


II. MODEL 

We consider the following Hamiltonian on the pyrochlore 
lattice: 

n^no + Ui (la) 

no = Yl JzzS/S^ - J± (5+.5- + S-S+) (lb) 
(V) 

Hi = ^ j±± (5+^+7y + 

(V) 

+ Jz±{S! {QjS+ + QS-) +^^J} (Ic) 

where Sf^ = Sf ± iS/ and S/, with /i = z, +, —, is defined 
in the local [111] coordinate frameattached to each of the 


We consider a s-MET treatment of the Hamiltonian of 
Eq. (1) which orders in a state with a fc = 0 ordering wave 
vector in a certain region of the parameter space centered 
around a dominant J± (Fs region). In this state, all the mag¬ 
netic moments in the system lie predominantly perpendicular 
to the local [111] direction and point in the same direction in 
the local coordinate system (see Fig. 1). These spin configu¬ 
rations define the Fs manifold.^® As a result, m-i = (Si), the 
on-site magnetization,^' and 0^ = tan“'(m//mf), the az¬ 
imuthal angle expressed in the local [111] coordinate system, 
are independent of the lattice site index i. We henceforth drop 
the index i of 0i and rrii for these fc = 0 spin configurations. 

As first noted in Ref. [16], in the Fs region, the s-MFT 
free-energy is independent of 0 and the system displays an 
accidental (7(1) symmetry within such a s-MFT description. 
However, since Ti. only has a global Zq symmetry, we expect 
that in a treatment of the problem that goes beyond s-MFT, 
this “artificial” (7(1) symmetry to be reduced to a Zq sym¬ 
metry in the paramagnetic phase which gets spontaneously 
broken in the ordered phase. In this context, we note that a 
high-temperature series expansion of the quantum model ( 1 a), 
with values {J±, J±±, Jz±, Jzz} appropriate for Er2Ti207,*® 
shows explicitly that such a Zq anisotropy develops in the 
paramagnetic phase upon approaching Tc from above. Also, 
a recent Monte Carlo study has investigated the critical behav¬ 
ior of a classical version of this model.We now present a 
Ginzburg-Landau (GL) symmetry analysis that allows one to 
anticipate the form of the lowest order fluctuation corrections 
to s-MFT that lift the (7(1) degeneracy. 
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A. Ginzburg-Landau (GL) symmetry analysis 

We start by defining the complex variable 

rrixy = exp(z(/)), (2) 

which corresponds to the magnitude and direction of the on¬ 
site magnetization m = (Si) in the local xy plane. 

For simplicity, we first consider a strictly local-xi/ state and 
assume = 0. The allowed terms in the Ginzburg-Landau 
(GL) free-energy can only be of even powers in m^y due to 
time-reversal symmetry, r. On the other hand, the rotation 
by 27r/3 about the cubic (111) axes, or C3 symmetry, forbids 
the existence of terms of the form rri^y -f unless n 

is a multiple of three. The resulting terms have the ability to 
lift the accidental U{1) degeneracy of the Ls manifold since 
they introduce a dependence of the GL free-energy on the az¬ 
imuthal angle (j), the orientation of m in the local xy plane 
(see Eq. (2)). Considering the effect of and r together, as¬ 
suming strictly xy order (m^ = 0), the lowest order term that 
breaks the U{1) symmetry must therefore have n = 6 and be 
of the form 


Ve{T)[mly +{mlyf], (3) 

where the anisotropy strength 775 depends on temperature, 
T. We note that such a sixth order term is dangerously ir¬ 
relevant for the 3-dimensional xy universality class.This 
means that while it does not affect the critical properties of the 
system,^^’^^ this term plays a crucial role in the selection of a 
specific long-range ordered state by lifting the U{1) s-MFT 
degeneracy at T < T^. Higher order terms ^ 

(n > 1) are also generated by fluctuations, but are even more 
irrelevant than /q at Tc- As an example, we discuss in Section 
IV the effects of competing /g and /12 terms at T < Tc as the 
amplitude \mxy\ grows upon cooling below Tc- 

If we now include the z component of the order parameter, 
i.e. niz 0, a C/(l) symmetry-breaking term arises at fourth 
order in the components of m in the GL free-energy. Such a 
term was only recently noted in a numerical study^' but whose 
microscopic and symmetry origins was not discussed. Again, 
based on the combined effect of the C3 and r symmetry op¬ 
erations, the degeneracy-lifting fourth order term in the GL 
free-energy has the form 

uj{T)mz [mly + [mlyf], (4) 

where, here again, the anisotropic coupling uj depends on T. 
Together, Eqs. (3, 4) identify the form of the two lowest order 
terms in the GL free-energy capable of lifting the degeneracy 
within the Eg manifold. 

It turns out that the sixth order term of Eq. (3) and the fourth 
order term of Eq. (4) have the same net effect in lifting the de¬ 
generacy beyond s-MET as can be shown by combining them 
into a single term in the GL free-energy, J^gl- Consider first 
the quadratic terms in of the form ro(m^ +'^y) +rim1. 
Here tq and ri are chosen to be different to emphasize the 
distinct criticality for the xy and z components of the order 


parameter m. Next, taking into account the terms of Eqs. (3, 
4), we have: 

-^GL = roiml + ml) + riml + + [mlyf] 

+r]&[{mly + {mlyf]+ F'^l[mxy,mz) 
TT^Ql{mxy,mz) ^ -. (5) 

This can be rewritten as 


2 I 2 N I / COs( 3 ^) . 2 


•T^GL = +mi) +ri{mz + 




2ri 


+ (2% - )|m 

2ri 


xy 


Tl 

® cos(6^) 


+J^'^l{mxy, mz) + J^Ql{mxy, mz) + ■ 


(6) 


where we used Eq. (2) to go from Eq. (5) to Eq. (6). In 
Eqs. (5,6), ^tid ^re fourth and sixth order terms and 
• • • represents terms that are of higher order in the compo¬ 
nents of m. We take ri > 0 to enforce no criticality for the 
mz component of m. Upon minimizing Eq. (6) with respect 
to mz, we obtain: 


=-a;|mxy|^cos(3()))/ri. (7) 


Hence, after having minimized with respect to m^, we 
have 


'5-^gl(()>) = 2?76 |TOxy|®cos( 6^), (8) 

with 

2j?6 = (2p6 - w^/2ri), (9) 

and where 5TGi^{(j)) is the anisotropic part of the C/(l) 
degeneracy-lifting term on which we focus. The minimum 
of bJ^GL is either cos{6(f>) = ±1 depending on the sign of 
fje. This happens for (j) = mr/3 or cj) = (2n + l)7r/6, with 
n = 0,1,..., 5 with these two sets of angles corresponding, 
respectively, to the -02 and the -03 states (see Fig. 1). So, when 
i?6 < 0 (76 > 0)- the minimum free-energy state is ip 2 (V'a)- 
The boundary between 7/>2 and 03 is thus determined by the 
real roots of pg- A similar discussion is invoked in Ref. [17] to 
describe the zero temperature 02-03 phase boundaries arising 
from quantum ObD. 

We note that the fluctuation corrections to s-MFT induce 
a local mz moment to the minimum free-energy state which 
was found to be strictly ordered in the local xy plane at the 
s-MFT level.^^ Since this moment is proportional to cos(30), 
only the 02 state can have a nonzero mz- This result is also 
expected based solely on symmetry considerations for the 02 
and 03 states (see Appendix B). The induced moment is, how¬ 
ever, small just below Tc since it is proportional to \mxy\^ 
and inversely proportional to ri which remains firmly posi¬ 
tive away from spontaneous Ising criticality at ri = 0. 

At the phenomenological GL level, Eq. (8) is the final result 
demonstrating how the C7(l) degeneracy at the s-MFT level 
may be lifted by the anisotropic terms of Eqs. (3,4). In this 
work, however, we are rather interested in exposing how the 
coefficient of the symmetry-breaking terms in the free-energy 
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of Eq. (6), fje, can be determined at a microscopic level when 
going beyond a s-MFT description. Specifically, we wish to 
explore the leading dependence of fyg upon the J±±, Jz± and 
Jzz anisotropic exchange couplings near Tc- 

III. METHODS AND RESULTS 

In Section III A, we use the E-TAP method to compute the 
phase boundary between the ip 2 and 1/^3 states determined by 
the function fje in Eq. (8) for T < . In this calculation, we 

consider for simplicity Jzz = 0 in Eq. (I).^^ In Section IIIB, 
we conduct a complementary numerical study of fluctuation 
corrections to the s-MET using cluster-MET (c-MET) which 
allows us to explore the effect of fluctuations at the level of 
one tetrahedron for both Jzz = 0 and Jzz 0. Since we are 
essentially interested in the E5 region where the dominant in¬ 
teraction is Jj- in Eq. (1), we express the rest of the couplings 
in units of J± and we henceforth denote the scaled (perturba¬ 
tive) interactions with lower case letters, i.e. j±± = J±±/J± 
and jz± = Jz±/J±- 

A. Extended-TAP Method (E-TAP) 

I. Method 

In a magnetic system with static magnetic moments, a given 
moment is subject to a local field due to its neighbors. In a 
s-MET treatment, the moment affects its own local field in¬ 
directly; this is an artifact of s-MFT. The so-called Onsager 
reaction field introduces a term in the s-MFT free-energy that 
cancels this unphysical effect to leading order. As was shown 
by Thouless, Anderson and Palmer (TAP),^® this reaction- 
field correction is particularly important in setting up a proper 
mean-field theory description of spin glasses with infinite- 
range interactions.^® Here, we present an extended version 
of the TAP method (E-TAP) first developed by Georges and 
Yedidia"*' for Ising spin glasses and which includes the lowest 
order fluctuation corrections originally calculated by TAP as 
well as those beyond. 

To proceed, we must first modify the E-TAP procedure of 
Ref. [41] to study the case of 3-component classical spins 
with anisotropic exchange interactions. In the E-TAP method, 
nonzero on-site fluctuations i.e. ^ are 

taken into account via a high-temperature expansion (small fj) 
of a Gibbs free-energy: 

G = In ^Tr[exp ( - /3H -f ^ A, • (5* - . 

( 10 ) 
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Here, rrii is the average magnetization at site i, rrii = (Si) 
and \i is a Lagrange multiplier which fixes rrii to its mean- 
field value. The high-temperature expansion introduces fluc¬ 
tuations about the s-MFT solution. Defining l3G{j3) = G{j3), 
the first two terms of the expansion in powers of /?, G(0)//3 
and G"(0), are the entropy and energy at the s-MFT level, 
respectively. The prime represents differentiation with re¬ 
spect to p. The higher order terms in the /? expansion of the 
Gibbs free-energy correspond to fluctuation corrections to the 
s-MFT free-energy that generate the terms that lift the t/(l) 
degeneracy. We aim to calculate the higher order terms (be¬ 
yond pP) in the P expansion of Eq. (10) that contribute to 
the degeneracy lifting term 2fjQ\mxy\^ cos(6(/)) of Eq. (8) to 
lowest order in the coupling constants j±± and jz± In other 
words, we are considering the selection near mean-field criti¬ 
cality and perturbative (in jz± and j±±) vicinity of the t/(l) 
symmetric portion of the theory (Hq in Eq. (la)). Since we 
are focusing on the expansion of G in Eq. (10), we write 
its correction 6G beyond the s-MFT solution as suggested by 
Eq.(8): 


SG = 2r]e{j±±,jz±)\mxy\^cos{6(j)). (11) 


As discussed in Section II A, fiQ{j±±, jz±) = 0 determines 
the phase boundaries between the 'ip 2 and the '03 states in the 
space of coupling constants j±± and jz±. We now proceed to 
calculate PeOiitj Jz±) with the E-TAP method. 

In the n**' order of the P expansion, factors of the form 
j±±fz± arise where k and I are positive integers and k+l = n. 
Each power of j±± and jz± contributes factors of and 
g±i0, respectively (see Eq. (Ic)), to the corresponding term 
in the P expansion. By a simple power counting of these fac¬ 
tors, one can pinpoint the terms that contribute to cos(60) and 
cos(30) that are necessary to compute due to their contribu¬ 
tion to uj and tjq and, consequently, to fje (see. Eq. (9)), and 
which arise at different orders in P in the E-TAP calculation. 
It is straightforward arithmetic to find what combinations of 
k and I in j±±ji± generate cos(6(/)) and cos(30) terms. We 
then calculate a; and pg in Efis. (5, 6) in terms of the micro¬ 
scopic couplings j±± and using the E-TAP method (see 
Appendix C). We obtain uj{j±±,jz±), Ve{j±±,jz±), and thus 
V6{j±±,jz±), all to lowest nontrivial order in j±± and jz±. 
These read as:®® 


^{j±±Jz±) = aoPj±±jz± + P‘^{aij'^^jz± + a2jl±) + /3^(a3j±±j^± -f a4j±±jf±) (12) 

and 

V6{j±±,jz±) = boP^j±± + + b2P^j±±jt± + b3P^jz±, 


( 13 ) 
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and, consequently, 

2%0±±, jz±) = Co/3^ j±± + + C 2 P'^j±±jt± + /3^(C3ji±jf± + Cij±±ji± + C5jf±)- (14) 


In Eqs. (12, 13, 14), a, (* = 0, • • • , 4), bj (j = 0, • • • , 3) and 
Cfe (fc = 0, • • • , 5) are numerical coefficients determined by 
the explicit E-TAP calculation too lengthly to reproduce here, 
but whose derivation is described in Appendix C. Considering 
the highest power of coupling constants j±± and jz± in Eq. 
(14), one needs to compute terms up to sixth order in the high- 
temperature expansion of G in Eq. (10). 


2. Results 

In order to obtain the numerical value of the a^, bj and Cfc 
coefficients in Eq. (14), we employ a diagrammatic technique 
to represent the terms in the /3 expansion and which consti¬ 
tutes the computational core of the E-TAP method. These di¬ 
agrams are composed of vertices and bonds.The vertices 
correspond to lattice sites covered by a diagram and the bond 
represents the interaction between the vertices that it connects. 
The details necessary to carry out the calculations using di¬ 
agrams are presented in Appendix C4. As mentioned ear¬ 
lier, for computational simplicity, we only consider the case 
of jzz = 0. Recalling Eq. (9), and noting that to the lowest 
order we can use the s-MET value of ri, ri = 2//3, we find: 


I 

2%(j±±, j.±) = 10-"^ X [-6.50/32^3^ + 57.3/33^ 26.5P^j±±jt± + 27.8P^jt± 

-7.95/33j4±j,\ - 4.60/34j3±j,\ + (15) 


do 

fli 

02 

03 

a4 

-0.593 

-1.13 

-0.451 

-0.938 

7.21 


TABLE I. Values of Oi in w(j±±, iz±) multiplied by 10^. 


bo 

bi 

^2 

63 

-0.322 

2.92 

1.42 

1.42 


TABLE 11. Values of bj in ??6(i±±,iz±) multiplied by 10^. 


In deriving Eq. (15), we used ui and tjq in Eqs. (12, 13), 
respectively, with the numerical values of at and bj coeffi¬ 
cients given in Tables I and II. Since fje{j±± ,jz±) = 0 deter¬ 
mines the boundaries between the tjj 2 and states, we need 
to compute the roots of Eq. (15). In order to be consistent with 
the above E-TAP derivation perturbative in jz± and j±±, we 
ought to only consider the lowest order term in the root of pg 
given by Eq. (15). We find that Eq. (15) has one real root 
which to, the lowest order in the coupling constants, reads: 

J±±~Poj1±3- (16) 

where the • • • represent higher order terms in jz± and po 
9.37/3 which determines the phase boundary between the ip 2 
and the states for jzz = 0, jz± -C 1 and /3 > /3c. The 
corresponding ip 2 -’ip 3 boundary for /3 = 0.253 > /3^^ = 1/4, 
providing an estimate of the boundary in the limit T T~, 
is shown Eig. 2 The “outer” s-MET boundaries of the overall 


r5 region for jzz = 0 is also shown in this figure. Outside 
the r5 boundary, the system orders in a Palmer-Chalker (PC) 
state'^46,60 Qj. ^ spjayed ferromagnet (SE)."^^ 

It was found in Ref. [17] that, at zero-temperature, quantum 
fluctuations yield three distinct phase boundaries for jzz = 0, 
a result confirmed in Ref. [32], even in the regime jz± < 1 
and j±± <C 1, for which the E-TAP results for the ordered 
state selection at T < Tc presented here apply. The combi¬ 
nation of the E-TAP results with those from Ref. [17] sug¬ 
gests, because of the different '02/'03 boundaries at T = T~ 
and T = 0, the possibility of multiple transitions between '02 
and 03 as the temperature is decreased well below Tc as was 
found for the model of Ref. [34]. Such a multiple-transition 
scenario constitutes an exotic variant of the more conventional 
ObD phenomenon since fluctuations select distinct long-range 
ordered states in different temperature regimes (T < Tc and 
T = O"*"). That being said, one should be reminded that 
since the E-TAP phase diagram of Pig. 2 was constructed 
on the basis of a lowest order fluctuation correction to s- 
MPT described by Eqs. (12,13,14) and for a classical ver¬ 
sion of T-L in Eq. (1), this discussion of multiple transitions 
is therefore only qualitative within the present E-TAP cal¬ 
culation. However, we expand further on this topic in Sec¬ 
tion IV within a Ginzburg-Landau theory framework. In the 
low-temperature regime, T <C Tc, a study that incorporates 
high-order magnon-magnon interaction at nonzero tempera¬ 
ture starting from the results of Ref. [17] would be of interest 
to explore this phenomenology further. This is, however, be¬ 
yond the scope of the present work. 
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Jz±/J± 

FIG. 2. (Color Online) The = 0 phase diagram of the model at 
T < Tc- The Ts region corresponds to the region enclosing the '02 
and 03 states taken together. The T 5 region is circumscribed by an 
outer (red) parabola j±± = — 2 and a horizontal red line at 

j±± = 2 obtained from the s-MFT calculation.*^ The phase bound¬ 
ary between the 02/03 states obtained from the E-TAP calculations 
(i.e. the roots of fj^ = 0 , see text) is represented by the black line. 
The phases outside of the boundaries of the Fs region are: a splayed 
ferromagnet (SF)"*^ canted from the [ 100 ] cubic direction and the so- 
called 04 or Palmer-Chalker (PC) state. Both the SF and the 
PC phases have k = 0 ordering wave vector. 


B. Cluster mean-field theory (c-MFT) 


It is interesting to ask to what extent a calculation, any cal¬ 
culation, that goes beyond s-MFT may reveal an ordered state 
selection at T < Tc- For this reason, we study in this sec¬ 
tion the effect of fluctuations numerically using c-MFT. This 
method incorporates fluctuations at the level of one tetrahe¬ 
dron by exactly diagonalizing the Hamiltonian corresponding 
to the tetrahedron and treating the interactions between dif¬ 
ferent tetrahedra at the s-MFT level. The c-MFT approach 
does not involve a perturbative scheme in powers of /3 com¬ 
pared with the E-TAP method where fluctuations are included 
through an expansion of the Gibbs free-energy in powers of 0 
(see Eq. (10)). In c-MET, fluctuations are limited to one tetra¬ 
hedron while, in the E-TAP method, fluctuations beyond one 
tetrahedron are included up to the range spanned by the dia¬ 
grams generating the sixth order terms in the high-temperature 
expansion of G (see Appendix C4 for details). One practi¬ 
cal advantage of the c-MET method compared to the E-TAP 
approach is that c-MET allows us to easily investigate the ef¬ 
fect of nonzero Jzz on the selection of the 02 versus the 03 
state. In what follows we first provide the details of the c-MET 
method and then present in Eig. 4 the results of the calcula¬ 
tions for various Jzz values. 


1. Method 


The c-MET that we use here may be viewed as a general, 
non-perturbative and rather system-independent way to obtain 
(numerical) results beyond s-MET. While the approach is not 
specific to the type of spins (classical or quantum) considered, 
we restrict ourselves to quantum spins 1/2 for a tetrahedron 
cluster. A pragmatic reason for using quantum spins here is 
to avoid the complicated eight integrals over the solid angle 
that correspond to the classical trace for classical spins taken 
as vectors of fixed length |5'i| = 1/2 and orientation defined 
by an azimuthal and a polar angle. To illustrate the method, 
we write the model Hamiltonian (1) in the compact form: 

E (17) 


with /i, lA = 2;, +, — 

and; 




J 22 

JzzkCij 

Jz±Qj 

'll 

-H 


-J± 



-J± 

J±±ilj 


Eor N sites on the pyrochlore lattice, there are N/A "up" and 
A^/4 "down" corner-sharing tetrahedra (see Eig. 3). The di¬ 
amond lattice, dual to the pyrochlore lattice, offers a simple 
representation of the tetrahedra that tessellate the lattice as el¬ 
ementary units. We note that the diamond lattice is bipartite (it 
is composed of 2 ECC sublattices, say A and B). As we con¬ 
sider ordered phases with a fc = 0 ordering wave vector, we 
assume in the c-MET method that all tetrahedra on sublattice 
A are equivalent and interact with each other only through the 
mean fields. We re-express model (17) in a more suggestive 
form as a sum over sublattice A (labelled with / or I') and 
take 0 j = 1, 2, 3,4 as the sublattice indices: 

nj _ \ ' \ ' cA* qv I \ ' \ ' qM qi/ 

(19) 


We then proceed to apply a s-MET approximation only on 
the second term of the Eq.(19): we neglect any fluctuations 
between the A tetrahedra while taking full account of fluc¬ 
tuations within the A tetrahedra. The self-consistently deter¬ 
mined mean fields introduced to decouple 

the inter-cluster bonds. Again, because we focus on the Fs 
manifold with fc = 0 ordering wave vector, we have transla¬ 
tional invariance of sublattice A which implies mJ'j, = mJ'j. 
By performing this approximation, Eq. (19) reduces to a sum 
over the A tetrahedra coupled together by the mean-fields: 


Hmf = 


E E 

I^A 


jllV 




A* 


( 20 ) 


Any thermodynamic average is readily computed from T^mf- 
In particular, it is straightforward to show that the physical 
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FIG. 3. The pyrochlore lattice, with spins located at the vertices of 
comer-sharing tetrahedra. There are two orientations of tetrahedra 
indicated by the black and red colored tetrahedra with the interaction 
paths represented by the colored lines. The centers of the tetrahe¬ 
dra form a diamond lattice (with sublattice A and B corresponding 
respectively to, say, the black and red tetrahedra). In c-MFT, the in¬ 
teraction in-between the black tetrahedra is treated at the mean-field 
theory level. 


ordered state built with such “ 2 -in/ 2 -out" spin configuration 
would have a nonzero net magnetization on each tetrahedron, 
as well as for the whole system since we consider a fe = 0 
propagation vector, which is forbidden for the Fs representa¬ 
tion under discussion. Therefore, the best that a positive jzz 
can do is to favor the state (uiz = 0 ) against t/i 2 (toz ^ 0 ) 
and the boundary shifts accordingly, as shown in Fig. 4. 




Jz±/j± Jz±/J± 


solution {mf} that minimizes the free-energy corresponds to 
rrpp = {SPp), where the thermodynamic average is taken with 
respect to the cluster mean-field Hamiltonian 'Hmf- While we 
focus on the = 1/2 case, we note that c-MFT could be 
used for | | > 1/2 and extended to clusters with more than 4 

spins. The method (using quantum spins) is only constrained 
by the computational limitations of the required exact full di- 
agonalization over the cluster considered.®' 


FIG. 4. (Color Online) The c-MFT phase diagram at T < Tl ob¬ 
tained for positive jzz = 0.5, 0.1, 0.01, 0. The region circumscribed 
by the red line is the 1 I 12 /'pz T 5 region. For a given phase boundary 
(fixed jzz) the area above (below) that phase boundary line corre¬ 
sponds to the 'i />2 (t/’s) phase. The phase boundaries for negative val¬ 
ues of jzz will be approximately a reflection of the presented phase 
boundaries with respect to the jzz = 0 line. Again, the phases out¬ 
side of the boundaries of the Ts region are: a splayed ferromagnet 
(SF) canted from the [100] cubic direction"'^ and the so-called 'ip 4 ~^ 
or Palmer-Chalker (PC) state.®® 


2. Results 

The critical temperature Tc is a function of the J±, Jz± 
and J±± exchange parameters. For every point in the phase 
diagrams (i.e. for a given set of exchange couplings jz± 
and j±±), we determine Tc by identifying the temperature at 
which there is a minimum of the free-energy and the develop¬ 
ment of a numerically non-zero on-site magnetization m. For 
definitiveness, we then take a temperature slightly below Tc, 
T = 0.9Tc, for each (jz±;/±±) determine which state 
{tp 2 or t/ 3 ) is selected. We record this state selection over a 
grid of points that span the F 5 manifold. 

The c-MFT tp 2 /tps phase boundary obtained for various jzz 
values is presented in Fig 4. The previous symmetry analysis 
in Section IIA of the F 5 manifold showed that a nonzero uiz 
component is only compatible with the ’ip 2 state. As a corol¬ 
lary and confirmation of this expectation, we observe that as 
jzz is varied, the 'ip 2 / 4’3 phase boundary position shifts. This 
can be understood physically in the following way. A neg¬ 
ative jzz favors a “4-in/4-out" hence non-zero and thus 
the tjj 2 state. Conversely, a positive jzz on a single tetrahe¬ 
dron alone would favor a “ 2 -in/ 2 -out" spin configuration and 
therefore disfavor a “4-in/4-out" configuration. However, an 


A comment is warranted regarding the flat (Jz±- 
independent) '!/’2/'03 phase boundary predicted by c-MFT for 
Jzz = 0 in Fig. 4. This is to be contrasted with the E- 
TAP results of Fig. 2 and the T = 0 quantum fluctuation 
calculations of Ref. [17] (see Fig. 3 therein), which both 
show a Jz±-dependent t/2-'03 phase boundary for Jzz = 0. 
This Jz± independent c-MFT result can be rationalized in the 
following way with the argument proceeding through a se¬ 
quence of steps. In the c-MFT approximation that considers a 
one-tetrahedron cluster, the ruz dependence of the free-energy 
comes only from the Jzz term because of the special role the 
C 3 symmetry plays when considering such a cluster. The rea¬ 
son is that every spin on a one-tetrahedron cluster is coupled 
to three idential “exterior” mean fields rrijj because of the 
fe = 0 tjj 2 or long-range ordered states considered. As 
a result, terms such as Qj'^j i) vanish due to 

the C 3 symmetry that imposes the necessary form for the Cij 
bond factors.'®’"'®’"'® The same argument of course applies for 
terms of the form Cij^j)- Consequently, the 

free-energy depends on explicitly only via the Jzz cou¬ 
pling since the combined dependence on Jz± and is elim¬ 
inated via the above argument. When Jzz = 0, all depen¬ 
dence of the c-MFT free-energy on disappears, and thus 
no-dependence on Jz± can remain, as found in Fig. 4. On 






















9 


the other hand, when J^z 7 ^ 0, the c-MFT free-energy de¬ 
pends on In that case, Jz± coupling the and com¬ 
ponents will, through the intra-tetrahedra fluctuations that 
are incorporated into the c-MFT calculation, renormalize the 
effective intra-tetrahedron zz component sublattice suscepti¬ 
bility and thus the net effect of inter-tetmhedm mean-field 
Jzzm^j j). As a result, for Jzz 7 ^ 0, a Jz± depen¬ 
dence of the phase boundary is observed for a single¬ 

tetrahedron c-MFT calculation (see Fig. 4). 

While the above argument for the Jz± independence of the 
boundary for the Jzz = 0 case applies for a single tetrahe¬ 
dron cluster c-MFT, we expect the boundary for larger clus¬ 
ters to depend on Jz±, even for Jzz = 0. In that case, Jz± 
may induce a dependence that would result in a behavior 
similar to E-TAP results for Jzz = 0. The reason is that, 
for larger clusters, the sites on the perimeter of the cluster 
may be coupled to the mean-field parameters rrijj coming 
from less than three sites related by C 3 symmetry, and which 
caused the inter-tetrahedron mean-field to vanish for a 4-site 
tetrahedron cluster as discussed above. Notwithstanding this 
caveat associated with the special symmetry of the 4-site tetra¬ 
hedron cluster, we nevertheless believe that the evolution of 
the boundary upon varying Jzz shown in Fig. 4 to be roughly 
qualitatively correct for Jzz 7 ^ 0 . 

IV. MULTIPLE TRANSITIONS AT T < Tc 

The results above, along with those of Refs. [17 and 32], 
suggest that in some circumstances the state selected at Tc 
may differ from the low-temperature phase selected by har¬ 
monic (classical or quantum) spin waves. For example. 
Ref. [34] found that, in a classical pyrochlore Heisenberg an- 
tiferromagnet with additional indirect Dzyaloshinskii-Moriya 
interaction, the state selected at T < Tj, is ^/'2 while the one 
selected at T = 0+ is V’ 3 - One might then ask what are the 
generic possibilities that may occur in xy pyrochlore magnets 
for which the anisotropic exchange terms position them in the 
degenerate Fs manifold. Considering the lowest order term in 
nixy, one has, as discussed above in Section IIA 

5J-gl(<^) = /6(T)cos(60). (21) 

fe{T) = 2fiQ\mxy\^ is a function of temperature, T, for fixed 
anisotropic exchange terms J±, Jzz, Jz± and J±±, with the 
sign of fe{T) dictating which state is selected at a given 
temperature. If feiT = O’*') stays of the same sign for all 
0 < T < Tc, only a single phase is realized below Tc- With 
solely this lowest order anisotropic term in SJ-qi,, the only 
other possibility is a first order transition between ip 2 and ip 3 
at some temperature T* < Tc when /^{T) changes sign at 
T = T*. However, as one gets far below Tc, higher order 
terms in \mxy\ can become of significant magnitude in J^gl- 
One may then extend the Ginzburg-Landau free-energy by in¬ 
corporating higher-order harmonics in the anisotropy poten¬ 
tial, as 

= Y. ‘^VeuiT)\mxyf^ cos(6n<^). (22) 

n=l 


Keeping only the two lowest order terms for illustration pur¬ 
poses, we have 

SJ^GLi^) = k{T) cos(6<^) + /i 2 (T) cos(12<^), (23) 

where /gn = ‘2,fjQn\'mxy\^'^■ For simplicity, consider first 
/12 = — 1 As /eir) varies from f&{T) < 0 to f%{T) > 0, 
a first order transition occurs when the minimum of SJ^g'l 
shifts discontinuously, in a first order transition, from cj) = 0 
(jp 2 State) to ^ = 7r/6 ('03 state) (see Fig. 5a). A more 
interesting behavior is in principle possible. Consider now 
/12 = +1. In that case, second order (Ising-type) transitions 
occur at /g = ±4|/i2|. For example, for /g < —4, the min¬ 
imum is at 0 = 0 (a ■02 state), and a second order transition 
to a /g-dependent angle 0 occurs at /g = —4. As /g contin¬ 
ues growing and become less negative (at fixed /12 = +1), the 
magnetic moment orientation 0 continues increasing untill an¬ 
other second order transition to the 03 state occurs at /g = 4-4 
(see Fig. 5b). While the direct 02 to 03 transition observed 
in the pyrochlore Heisenberg antiferromagnet with indirect 
Dzyaloshinskii-Moriya interaction^'* suggests that this system 
belongs to the first case, the second scenario is perhaps not 
excluded in more complex models in which both temperature 
and applied magnetic field are varied simultaneously, or when 
the magnetic Er^+ sites are diluted by non-magnetic ions.®^’®^ 



FIG. 5. (Color Online) (a) Sequence of (5 J^gl( 0) with /12 = —1 
and 06= -2, -1, 0, 1 and 2 (from top to bottom). There is a first-order 
transition at /g = 0. (b) Sequence of 57 ^gl( 0) with /12 = 1 and 
/6=-6, -4, -2, 0, 2 4, and 6 (from top to bottom). In mean-field theory, 
there are Ising-type second-order transitions at /g = ±4|/i21. 
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In the above discussion, we have implicitly assumed that 
the anisotropic terms [mz{mxy)^ + c.c.]^ and [{rrixy)^ + 
c.c.]^ which give terms of the form \mxy\^^ cos{6n(j>) in the 
Ginzburg-Landau theory once rriz has been eliminated (see 
discussion in Section II A), are generated by thermal fluctu¬ 
ations beyond the Ginzburg-Landau free-energy derived by 
a s-MFT treatment of the microscopic model TL in Eq. (1). 
However, it is in principle possible that virtual crystal field 
excitations mediated by the bare multipolar 

interactions between the rare-earth ions would generate^ 
an effective pseudospin-1/2 Hamiltonian more complex than 
the one given by Eqs. (la), (lb) and (Ic) and involve mul¬ 
tispin (e.g. ring-exchange like) interactions capable of lift¬ 
ing degeneracy at the classical level without invoking order- 
by-disorder.However, one naively expects those interac¬ 
tions generated by VCEE to be signigicantly smaller com¬ 
pared with the J±, J±±, Jz± and Jzz of Eq. (la), as stated in 
Ref. 16. Whether those terms are truly inefficient in compet¬ 
ing with order-by-disorder, once the commonly large prefac¬ 
tors of combinatoric origin arising in high-order perturbation 
theory have been taken into account, must await a detailed 
calculation. 


V. DISCUSSION 

In this work, we first used an extended TAP (E-TAP) 
method to analytically study the problem of order-by-disorder 
(ObD) near the critical temperature of a general three- 
dimensional xy antiferromagnetic model on the pyrochlore 
lattice with the Hamiltonian of Eq. (1) with Jzz = 0. We 
focused on the Es manifold which is 17(1) degenerate at the 
standard mean-field theory (s-MET) level. The fluctuations 
corrections to the free-energy beyond s-MFT were organized 
as an expansion in powers of the inverse temperature, /3, and 
to lowest orders in the coupling constants j±± = J±±/J±, 
jz± = Jz.±l J± and the on-site magnetization m. We estab¬ 
lished that in different parts of the E 5 region, the ■02 and 
states are the only minima of the free-energy selected by fluc¬ 
tuations up to the lowest order t/(l) symmetry-breaking terms 
considered. The phase boundary between ^2 and ips can then 
be obtained by finding the real roots of Eq. (14) in terms of 
J±±/J± and Jz±/J±- We also numerically studied the ObD 
mechanism in the 3D-xy pyrochlore system a using cluster 
mean-field theory (c-MFT). Using this method, we obtained a 
phase boundary between ^/>2 and '03 states for a variety of Jzz 
values shown in Fig. 4. 

Using a Ginzburg-Landau theory, along with E-TAP and c- 
MFT calculations, we predict that for a state ordered in the 
local xy plane at the s-MFT level, fluctuations can induce a 
small out-of-xy-plane niz component of the on-site magne¬ 
tization. This fluctuation-induced local z component of the 
magnetization is only compatible with the ip 2 state. We expect 
the size of this moment to be small, since it is proportional to 
p <C 1 for T < (see Eq. (7)). Yet, on the basis of 
the c-MFT results, we might anticipate that nonzero niz has 
an important effect on the phase boundary between the 1/2 and 
the '03 states for Jzz 0, as illustrated in Fig. 4 (see also the 


caption of the figure). 

By considering the single phase boundary between the 02 
and 03 states for T < along with the multiple 02-03 
boundaries found at T = O'*" in Ref. [17], for j±± <C 1 
and jz± 1, the regime for which the perturbative E-TAP 
solutions above apply, one may expect to find multiple phase 
transitions between 02 and 03 states upon decreasing tem¬ 
perature from T < Tc to T ~ O’*'. While possible in prin¬ 
ciple and found in a numerical study,® such multiple transi¬ 
tions in real xy pyrochlore magnetic materials have not yet 
been reported. Since in the E-TAP calculations the j±± and 
jz± couplings were treated perturbatively, the E-TAP phase 
boundary in Eq. (16) is valid for small j±± and jz± and 
would thus need to be modified for non-perturbative j±± and 
jz±, that is, closer to the boundaries of the Es region with 
the splayed-ferromagnet (SE) phase (see Eig. (2)). A case in 
point is the Heisenberg pyrochlore antiferromagnet with indi¬ 
rect Dzyaloshinskii-Moriya interaction.®^ This model displays 
a transition at r+ from a paramagnetic state to 02 , followed at 
T~ < by a transition from 02 to 03,®® and for which the 
corresponding anisotropic exchange J±, J±±, Jz± and Jzz 
are such that this model lives right on the Es to SE classical 
phase boundary.*^ 

The difference between the phase boundaries obtained by 
E-TAP and c-MFT for Jzz = 0 (see Figs. 2 and 4) arises 
from the different range of fluctuations considered in these 
two methods. In c-MFT, because of the exact diagonaliza- 
tion of the Hamiltonian on a single tetrahedron and the s- 
MFT treatment of the inter-tetrahedra couplings, the fluctu¬ 
ation corrections considered are short-ranged (limited to the 
nearest-neighbors). However in the E-TAP method, by con¬ 
sidering the higher order terms in the /3 expansion, specifically 
/3^ and beyond where larger (more extended) diagrams appear 
in the calculation (see Appendix C, Section C 4), fluctuations 
beyond nearest-neighbors are included. It would be of interest 
to benchmark these arguments by performing c-MFT calcula¬ 
tions for larger clusters. 

The E-TAP corrections to the s-MFT free-energy could also 
be computed for the Jzz ^ 0 case. The procedure for consid¬ 
ering Jzz 0 0 is conceptually the same as for the Jzz = 0 
case for which, using the E-TAP method, we calculated the 
terms in the high temperature expansion of Eq. (10) that con¬ 
tribute to the degeneracy-lifting cos(30) and cos(60) terms in 
the free-energy. To obtain these terms to the lowest order in 
the coupling constants J±±, Jz± and Jzz, one would need to 
consider the high temperature expansion of G up to the sev¬ 
enth order in 0 and calculate the degeneracy lifting terms fol¬ 
lowing the power-counting prescription explained above Eq. 
(12). However, the number of terms of the form j\,zj±±jz±, 
with fixed I + m + n value I, m and n positive integers, and 
thus the number of diagrams to be computed, proliferate sig¬ 
nificantly as to make this calculation a significant undertaking. 
The Jzz ^ 0 case thus stands on its own as an independent fu¬ 
ture study. 

Finally, we remark that variants of the E-TAP method pre¬ 
sented in this work could be applied to models other than 
the pyrochlore structure to analytically investigate the role 
of fluctuations beyond s-MFT in selecting a specific long- 
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range ordered state close to the critical temperature. One ex¬ 
ample is the case of Heisenberg spins on face centered cu¬ 
bic (FCC) lattice interacting via magnetostatic dipole-dipole 
interaction.®®’®^ Generally speaking, one might expect that a 
consideration of an E-TAP analytical description of fluctua¬ 
tions at the microscopic level may help shed light on the role 
of individual symmetry-allowed interactions in the ordered 
state selection near the critical temperature over a broad range 
of frustrated spins models for which an accidental degeneracy 
exists in a standard mean-field theory description. 
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Appendix A: 1 P 2 and ip 3 states 


state), the improper rotation by tt/2 (S 4 ) about the same cu¬ 
bic axis, and two plane reflections cri and (T 2 with respect to 
planes spanned by the cubic axis and one of the two tetrahe¬ 
dron bonds perpendicular to the cubic axis. The generators of 
the little group of '02 states include TS 4 , rai and tct 2 where 
T is the time-reversal symmetry operation. 

Considering the action of the symmetry operations in the 
02 ’s and ' 03 ’s little groups on the possible configurations with 
a finite onsite rriz moments on a single tetrahedron (i.e. all- 
in/all-out, 2-in/2-out, 3-in/l-outand l-in/3-out where “in” and 
“out” indicates whether rrii^z is positive or negative on site i), 
only all-in-all-out configuration is invariant under the symme¬ 
try operations of the ■02 s little group while none of the above 
configurations are invariant under the 03 ’s little group sym¬ 
metry operations. As a result, only the 02 state can possess a 
finite onsite ttIz moment which displays an all-in-all-out con¬ 
figuration on a single tetrahedron and therefore does not pro¬ 
duce a net magnetic moment on a tetrahedron. 


Appendix C: Extended TAP Method (E-TAP) 

In this Appendix, we first derive the general form of the 
E-TAP corrections up to sixth order in /3. Next, we illustrate 
our method for calculating the terms in the (3 expansion of the 
Gibbs free-energy, G, of Eq. (10) by focusing on two calcu¬ 
lation examples. Einally, we discuss the general properties of 
the diagrammatic approach employed to carry out the compu¬ 
tation of the various terms entering the E-TAP calculation. 


The 02 states have spin configurations with a A: = 0 or¬ 
dering wave vector and with the following spin orientations 
on a tetrahedron expressed in the global (Cartesian) reference 
frame 



(Ala) 


(Alb) 


(Ale) 


(Aid) 


Here the subscripts correspond to the four sublattice labels in 
the pyrochlore lattice and the superscripts refer to different 
symmetry-related 02 states, {ef} and {e?}, i = 0, • • • ,3, can 
be obtained from Eq. (Al) by C 3 (tt/S) rotations with respect 
to the ( 111 ) directions. 03 states can be obtained from the 02 
ones, by a 7r/6 rotation of each spin about its local [111] axis. 


Appendix B: Symmetry groups of 02 and 03 states 


1. Derivation 

The E-TAP method is based on a perturbative expansion 
of the Gibbs free-energy as a function of inverse tempera¬ 
ture, P, beyond that given by s-MET. There are several ways 
of deriving these corrections.®* Here, we employ and extend 
the method presented by Georges and Yedidia^' for Ising spin 
glass which provides the corrections due to fluctuations about 
the s-MET solution, order by order in 0. In this part of the Ap¬ 
pendix, we focus on a temperature regime close to the mean- 
field transition temperature, T < yMP, derive the E- 
TAP corrections for classical Heisenberg spins of fixed length 
| 5 | = 1 / 2 . 

The Hamiltonian in Eq. (1), TL, can be written as: 

H = (Cl) 

i<j 

Greek labels represent Cartesian coordinates and implicit 
summation over repeated superscripts is used. A Taylor se¬ 
ries expansion of Eq. (10) in powers of /? reads: 


The symmetries of either 0^ state (i = 2, 3) form a group 
known as the little group of the corresponding state. The gen¬ 
erators of the 03 little group include the C 2 rotation by tt about 
one of the cubic x,y or z axes (depending on the particular 03 


G(/3) 




dGjP) 

dp 


1 d^G(p) 




(C2) 
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where G(/3) = We define. 


(C3) 

ij 


where 

6SJ = - mj" (C4) 

are spin fluctuations about the s-MFT solution. As shown in 
Ref. [41], the derivatives of G{P) with respect to /3 can be 
evaluated in terms of expectation values of powers of U and 
Tn, that read; 


dmm 

dl3 

d^iPGm 

a /32 

d^^Gm 

d\PGm 

d/3* 

dHm/3)) 

d/3^ 

d^mm 

5/36 


(n), 

(U^), 

+3{U'^f -3{U^T2), 

(U^) - 1Q{U‘^){U^) - 3{U'^T3) + 7{U^T2) + 6(C/T|), 

-(C/®) + 15{U^){U^) + 10{UY - S0{U^f - 12{U'^T2) 

+10{U^T3) - S{U^T4) - 27{U^T^) + 18{UT2T3) 

-6{UT3){U^) + 51{U^){U^T2) + 6 (C/ 2 )(r|) - 6 (T|), 


(C5a) 

(C5b) 

(C5c) 

(C5d) 

(C5e) 

(C5f) 


and where T„ is defined as 

(C6) 

i 

Since dX/d/3\i3=o = hi*^ where fif = J2j,u 
terms involving involving T„ come from considering fluctua¬ 
tions of the local mean-field. On the other hand, the C/" terms 
take into accont the fluctuations of the on-site magnetization 
(see Eqs. (C3, C4)) within the ensemble set by the s-MFT 
solution. The (• • •) above denotes a thermal average. For a 
general observable O, (O) is given by; 

^ Tr[0 exp ^ + " "^»))] .(. 7 . 

rr[exp ( - I3H + Xi ■ (Si - m*))] 


is nonzero only if there is no site label that appears only once. 
For example, averages of the following form have a nonzero 
contribution; 


{SSYSSYSSYSS°^^SSY) = {SSYSSY){SSYSSYSSY) 

(C9) 

The expectation values above can be calculated using the 
self-consistent s-MFT equations for 3-component classical 
spins which are given by the Fangevin function 



coth(|Aj|) 



(CIO) 


Consequently 


{ssrssf) 


dmf _ ap 

aAf 


(Cll) 


The first two terms in Fq. (C2) correspond to the s-MFT free- 
energy while the higher order terms in /3 provide corrections 
beyond s-MFT. Calculating the expectation value of powers of 
[7 at /3 = 0 reduces to the evaluation of mean-field averages 
of the form; 


(SS^SS^SS^) = (C12) 

and in general 


{SS^^^SSZ^---SSI-)mf, (C 8 ) 

where *„ represents the site label and an represents a Carte¬ 
sian coordinate, n is the number of 6S factors in Fq. (C 8 ). 
From now on, we drop the MF subscript for compactness of 
notation. For n = 1, the expectation value in Fq. (C 8 ) is zero 
due to the relation rrii = (Si). For n > 2, however, Fq. (C 8 ) 


{SS^SSY ■ ■ ■ SSY) = 


d/SS^SSY ■■■SSY~") 


. (C13) 


Since we are interested in a temperature range close to Tc, 
Fq. (CIO) can be expanded for small jA^j; 


, Af (|A.|)2 2(|A.|)^ 

* 3 ^ 15 315 


(C14) 


15 


315 
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The expansions of Eqs. (C11,C12, C13) at T < can be 
computed by differentiating Eq. (C14) with respect to differ¬ 
ent components of the vector Lagrange multipler. 

In the rest of this section, we focus on the specific prob¬ 
lem of the s-MET U{1) degeneracy of the Es manifold dis¬ 
played by the Hamiltonian of Eq. (1) in a wide range of 
the {J±, J±±, Jz±, Jzz] anisotropic exchange parameters.'^ 
Below, we calculate the degeneracy-lifting contributions from 
(C/^) and (t/^T2) in Eqs. (C5b, C5d) to demonstrate the gen¬ 
eral idea of the method. Higher order terms can be calculated 
using a similar procedure. 


2. ([/^) Calculation 


Using Eq. (C3), we have 


(£/") 


bi ji> (*2^2) 


(C15) 

where the summations are performed over lattice sites and 
summation is implied for Greek superscripts. To proceed, we 
need to specify all the possible pairings of the factors 5Si^, 
SS,, (f = 1,2) in the 

Considering the description provided above Eq. (C9), the 
only nonzero pairings of site indices in Eq. (C15) are ii = 12 , 
ji = j 2 and ii = j2, ji = 12 , with the constraint ii 7^ ji and 
*2 ^ J2 imposed by the Hamiltonian; = 0 for all a and 
/3. This leads to: 


(C/2) = ^ X! (CI 6 ) 

(*b) 


The computational complexity of the method requires the 
usage of a diagrammatic approach, which we now introduce 
by considering, for example, the calculation of Eq. (C16). 
In this equation, for a given i, j, at and j3t with t = 1,2, 
{6S'^'^SS^^){SSj^SSj^) can be represented by a 
diagram of the form illustrated in Eig. 6a. In this figure, the 
vertex labels match the labels of the lattice sites that the dia¬ 
grams cover. Each vertex represents an average of the form 
((55“^ • • • (55“*) where t is the number of bonds (solid lines) 
connected to that vertex. In the case of Eq. (Cl6) as indicated 
in Eig. 6a, the averages are (b5“*(55“^) and {SS^^SS^^). The 

bonds represent the elements of the coupling matrix, It 
is straightforward combinatorics to take into account only the 
terms that contribute to the degeneracy-lifting factors, cos(3()>) 
or cos(6())). We next proceed to demonstrate this point. 

Based on the symmetry analysis of Section II A, we note 
that since ([/^) involves four powers of b5“*. (C/^) can 

only contribute degeneracy-lifting terms of the general form 
(X mz\mxy\^ cos{34>) and where the proportionality factor 
is generated by the anisotropic couplings j±± and jz± in 
Eq. (Ic). We recall that j±± and jz± each contributes a 
factor and in the power counting method. The 

—i(j) and —2i(j) corresponds to presence of SS~ and {SS~)‘^ 
in Eq. (C16), respectively, which in turn, implies the pres¬ 
ence of J-““ = j±±^*j or J^~ = jz±Qj matrix elements 



FIG. 6. a) Diagrammatic representation of the term 

m-Eq. (C16). b) All di¬ 
agrams corresponding to at, Pi = z,+, i = 1,2 with 2 : appearing 
only once among the Greek superscripts. The same number of 
diagrams are present for the case where ai,Pi = z,—. Again 2 
appears only once. 


in ). So for a given i 

and j, only the following combination of terms can generate 
cos(3<))): 


J++J^;{SSt6St){SS]SSp + 
J~J^-{SS-SS-){SS]SS-), (C17) 

where the first term generates the contribution to cos(3<))) 
while the second term generates All possible ways of 

generating the factor are illustrated in Eig. 6b, which is 
the same as the number of ways of generating Ulti¬ 

mately, Eq. (C16) can be rewritten as: 

(t/') = 2J±±Jz±[{SS^6S+){SS+SS+)J2l^JC^J+i^■c)■ 

u' 

(CIS) 

Based on Eqs. (Cll, C14) and recalling from Section II that 
m and (j) do not require site indices in the Es manifold, we 
have dropped the site labels of 6S’s in Eq. (C16) when writing 
Eq. (CIS). The lattice sum, can be carried out 

using a computer program for different lattice sizes with linear 
dimension L. Up to sixth order in the [3 expansion, the lattice 
sums per site for different terms in Eqs. (C5) are independent 
of L for L > 2 . Accordingly, we have 

^ (C19) 

ij 

Using Eqs. (Cll, C14 ), the averages in Eq. (CIS) can be 
written as 


{SS^SS~^) = —mzm+ + ■ ■ ■ (C20) 

1 8 

{SS+SS+) = + (C21) 
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FIG. 7. For details, see the text in Appendix C4. We also note that 
the contribution of disconnected diagrams (e.g. the diagram at the 
bottom right comer) in Eq. (C5) adds up to zero. 


We note that Eq. (C22) contributes to uj in Eq. (6) which, in 
turn, is necessary to evaluate rye(j±±, jz±) in Eq. (11). All 
other terms in Eq. (C5) involving solely powers of U and 
no Tn terms can be calculated in a similar way. The number 
of diagrams increases as one considers higher order terms in 
the /3 expansion of Eq. (10). As a result, finding the number 
and type of nonzero average of the form of Eq. (C8) is most 
easily done using a computer program. The details of this 
type of calculations are presented in Appendix C 4. Some of 
the diagrams that appear at higher order in the /3 expansion 
are illustrated in Eig. 7. 


3. {U^T 2 ) Calculation 


Due to the presence of T„ in Eq. (C6), in this case n = 2, 
averages that involve T„ need to be carried out slightly dif¬ 
ferently in comparison with averages containing only powers 
of U. The difference comes from the presence of the factor 
d^^Xi/dl3^ in Eq. (C6). Considering Eq. (C5a) and the rela¬ 
tion dX/d(3\i3=o = we obtain the following relation 


where to the lowest order of interest in m which, in this case, 
is the fourth order, we kept A“ ~ 3m“ and neglect all higher 
order terms. Einally Eq. (C18) gives 


{U'^)/N = (-0.96mzTO+ -F h.c.)j±±jz± 

= -^■Q‘^j±±jz±mz\m^y\^cos{3(j)). (C22) 

I 


_ d{H) _ d'^{/3G{/3)) 

dp dm°‘ dm'^dp ’ 

Erom Eq. (C23), one can write 

i9"A“ _ d^+\pG{p)) 

9/3” 9m“9/3” 


(C23) 


(C24) 


Using Eq. (C6), {U‘^T 2 ) can therefore be written as: 


{U^T2) 


1 _ _ _ 

SS'^^). 

ilh izjz k 


(C25) 


Using Eq. (C24), we have 


a^Af 93(/3G(/3)) 9(C/2) 


9^2 


9to“ 9/32 9m“3 


(C26) 


The term in Eq. (C25) can be 

dealt with as described previously in Appendix C 1 . The out¬ 
come reads 

{6S^p <55'“/ ) = 4(b5“/ b5“/) (<55/^ <55f ^ <55^^), 

(C27) 

where the factor of 4 comes from the number of different ways 
of pairing site indices yielding non-vanishing results. On the 
other hand, 

(C/2) = J“/'^V“='^^(<55r<55r)((55f ,55f), (C28) 


and in turn 


dm ^^ V m 

a™ “3 0™0!3 0\a3 '“‘^.7 “‘^.7 9 


9m“^ 9m“^ ^ */ 


dXf 


where 


dm. 


(C29) 

~ 3. Here, we ignored higher order terms 
which do not contribute to the degeneracy lifting terms of 
fourth or sixth order in components of m in Eq. (C25). The 

^ expression can be calculated using Eq. (Cl2). 
Substituting Eqs. (C27, C29) in Eq. (C25), we obtain the fi¬ 
nal expression which can be represented by a “fused” diagram 
shown in Eig. 8. This diagram has a new type of vertex repre¬ 
sented with a red square. This vertex is labeled with only one 
site index and its mathematical expression is: 


(<55/M5f ^ <55f ^) (<55/^ <55^^ 5Sl ^), 


(C30) 





























15 


where again there is a sum over repeated Greek superscripts. 
As a result, the final expression reads: 


J 


(C/^Ta) 




Y. SS^^){6S^^6S^^SS^^){6S^^){6Sf^6Sf^). 


(C31) 


From this point on, one can adopt the procedure presented pre¬ 
viously for the ([/”) terms to obtain the final result. All other 
averages that contain r„ e.g. and 

where tii, rrii with i = 1,2 are natural numbers, can be calcu¬ 
lated following a similar procedure. 



FIG. 8. “Fused” diagram corresponding to {U^T 2 ). The 
vertex represented by the (red) square represents the following: 

{SS^^ >( 55 ^ 3S2^ SS2^). 


4 . On Diagrams 

In this subsection, we make a few comments about the dif¬ 
ferent diagrams that appear in the E-TAP expansion. 

First we focus on diagrams corresponding to ([/") with 
n = 2, 3, • • •, in Eq. (C5). The diagrams for n = 3,4 are 
illustrated in Fig. 7. In this figure, blue circles represent a 
given lattice site and the number of lines (bonds) connected 
to the circles is the number of paired SS’s at that lattice site. 
The degree of the vertices is indicated above each diagram 
and is written in the form of (ai, • • • , a^) where m is the 
number of vertices and is the number of bonds connected 
to vertex i. The number of times that each diagram appears 
in the process of pairing SS’s, referred to as diagram count, 
is indicated in parentheses to the right of each diagram in 


Fig. 7. The diagrams are distinguished by the degree of their 
vertices and their adjacency matrix eigenvalues.^® The adja¬ 
cency matrix M is a m x m matrix with matrix elements 
Mij (0 < i,j < m) corresponding to the number of lines 
connecting the f-th vertex to the j-th vertex. We enumerate 
these diagrams using a computer program in which we gen¬ 
erate all possible diagrams given a certain number of vertices 
of a given degree. Then, we remove all diagrams that include 
any onsite interactions which are forbidden since there is no 
onsite interaction in the Hamiltonian of Eq. (1). For the re¬ 
maining diagrams, we build their adjacency matrix which we 
then diagonalize. M is a symmetric matrix and thus has a set 
of real eigenvalues which defines a graph equivalence class. 
Graph with the same adjacency matrix eigenvalues correspond 
to graphs with the same topology, i.e. they are isomorphic.®® 
The number of graphs with a given topology is noted in paren¬ 
theses to the right of the graph in Fig. 7. 

The fifth and sixth order terms of the E-TAP expansion in /3, 
contain, respectively, 9 and 26 types of connected diagrams. 
We note that the contribution of disconnected diagrams in Eq. 
(C5) adds up to zero.®® 

The so-called fused diagrams appear in terms of the form 
(C/"*T„) and (TmT„), where m,n > 2 are positive integers. 
They are constructed by fusing the diagrams similar to the 
ones in Fig. 7 together. An example of such a fused diagram 
and the corresponding details of its definition is given in the 
caption of Fig 8. 

We note that some of the diagrams do not cover lattice sites 
beyond one tetrahedron, for example the triangular diagram 
and the two site diagram in the top row of Fig. 7. Only, di¬ 
agrams of this type contribute to fluctuations incorporated in 
the cluster mean-field theory (c-MFT) calculations. However, 
their effect is incorporated to all orders in /3 in the c-MFT cal¬ 
culation. 
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